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A MODEL FOR RATE-DEPENDENT HYSTERESIS IN PIEZOCERAMIC MATERIALS 

OPERATING AT LOW FREQUENCIES 

RALPH C. SMITH*, ZOUBEIDA OUNAIES 1- , AND ROBERT WIEMAN* 


Abstract. This paper addresses the modeling of certain rate- dependent mechanisms which contribute 
to hysteresis inherent to piezoelectric materials operating at low frequencies. While quasistatic models 
are suitable for initial material characterization in some applications, the reduction in coercive field and 
polarization values which occur as frequencies increase must be accommodated to achieve the full capabilities 
of the materials. The model employed here quantifies the hysteresis in two steps. In the first, anhysteretic 
polarization switching is modeled through the application of Boltzmann principles to balance the electrostatic 
and thermal energy. Hysteresis is then incorporated through the quantification of energy required to translate 
and bend domain walls pinned at inclusions inherent to the materials. The performance of the model is 
illustrated through a fit to low frequency data (0.1 Hz - 1 Hz) from a PZT5A wafer. 

Key words, rate-dependent hysteresis model, piezoceramic materials 

Subject classification. Structures and Materials 

1. Introduction. The majority of currently employed models for hysteresis in ferroelectric materials 
are based on the assumption of static or quasistatic operating conditions. However, it has long been rec- 
ognized that the polarization which is generated at a given field strength is dependent upon the rate at 
which the field is cycled. Hence while quasistatic models may be suitable for initial material characteri- 
zation in certain applications, the incorporation of rate-dependence in the models is necessary to quantify 
the material behavior through its full operational range. In this paper we characterize the low frequency 
rate-dependent hysteresis in piezoceramic materials through a model comprised of two components: (i) A 
frequency-dependent anhysteretic model developed from Boltzmann principles, and (ii) Algebraic and ODE 
relations which quantify reversible and irreversible polarization changes due to the bending and translation 
of domain walls. 

The dependence of the polarization, in barium titanate, on the rate at which the field is cycled was 
discussed in detail by Landauer et al. [11] with significant reference to observations made by Merz [12]. 
These references illustrate that the coercive field at which the polarization switches direction is dependent 
upon the time allowed for switching. As illustrated in Figure 2 of [11], the rate at which the field cycles also 
significantly affects the shape of the polarization curve, with the maximum polarization attained at a fixed 
field strength decreasing with increasing frequency. 

To illustrate these phenomena in the context of piezo ceramic materials, data collected from a PZT5A 
wafer for input fields ranging from 0.1 Hz to 1 Hz is plotted in Figure 1.1. A comparison between the 0.1 Hz 
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Fig. 1.1. Rate- dependent hysteresis measured in a PZT5A wafer for input fields cycled from 0.1 Hz to 1.0 Hz. 

quasistatic data and the 1 Hz data illustrates a decrease in both the coercive field and the polarization as 
the frequency increases. This indicates the necessity of considering the effects of rate-dependence, even at 
very low frequencies. 

Numerous modeling strategies have been employed to quantify hysteresis in piezoelectric and ferroelectric 
materials. These include microscopic theories applied at the lattice or grain level [13], macroscopic theories 
based on empirical observations [4, 6, 7, 8, 9], and semi-macroscopic models which combine energy relations 
with macroscopic averages to quantify the bulk behavior of the material [3, 10, 15, 16, 17]. While some of 
these models incorporate frequency-dependence (e.g., [13, 17]), the majority of current analysis is directed 
toward static or quasistatic hysteresis phenomena. Moreover, a basic tenet underlying the construction of 
Preisach models is the assumption that the hysteresis behavior is rate independent [2], Hence this approach 
will not accommodate variable frequencies of the type illustrated in Figure 1.1 using a single parameter set. 

The model developed here is based on the approach employed in [14, 15, 16] for quasistatic regimes. The 
anhysteretic polarization is modeled first through the balance of electrostatic and thermal energies using 
Boltzmann principles. In the original formulation of the model, the anhysteretic polarization was formulated 
under equilibrium conditions between the field and polarization. Hysteresis was then incorporated by quanti- 
fying the irreversible and reversible changes in polarization due to domain wall movement. The modification 
of the model considered here is based on probabilistic arguments which ascertain the rate-dependence of the 
anhysteretic polarization. The resulting model yields the decrease in polarization observed in piezoceramic 
materials as the drive frequency increases. Hence it quantifies the effects observed in Figure 1.1. Finally, 
this model reduces to the quasistatic models in [14, 15, 16] when frequencies are limited to zero. 

The quasistatic model from [14, 15, 16] is summarized in Section 2 to illustrate the methodology and 
to indicate necessary modifications. The rate-dependent model is then developed in Section 3 and the 
performance of the model is illustrated in Section 4 through a comparison of the model prediction with the 
experimental data plotted in Figure 1.1. 

2. Quasistatic Model. To illustrate the modeling approach and indicate components which must be 
modified, we summarize the model developed in [14, 15, 16] for hysteresis in quasistatic regimes. As indicated 
previously, the model is comprised of two components. The first models the anhysteretic polarization which 
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is due to polar switching in addition to domain rotation at high field levels. Under certain conditions, the 
anhysteretic polarization is multivalued and hence incorporates a form of hysteresis. The transition between 
remanence and the coercive point is typically steeper than that observed in most ferroelectric materials, 
however, due to their polycrystalline nature and the inhibition of domain wall movement due to inclusions 
inherent to the materials. The latter effects are quantified through the consideration of energy required to 
bend and translate domain walls pinned at inclusions in the material. The combined model characterizes the 
nonlinear and hysteretic relation between the input field E and the polarization P generated in the material. 

The anhysteretic polarization for both quasistatic and dynamic regimes is modeled through the balance 
of the electrostatic and thermal energy using Boltzmann principles. For a dipole p in an electric field E, the 
potential energy is 

(2.1) S = -p • E = —pE cos(0) 

where p = |p| and E = |E|. The thermal energy at temperature T is given by £t = k B T where k B denotes 
Boltzmann’s constant. The probability that a dipole occupies the energy state S is then specified through 
Boltzmann statistics as 

( 2 . 2 ) n(£) = Ce- £ t kBT 

where the parameter C is chosen to ensure that integration over all possible dipole configurations yields the 
total number N of moments per unit volume. The assumptions specifying possible moment orientations 
determines the form of the anhysteretic model. 

The simplest model results from the assumption that dipoles can be oriented only in the direction of the 
electric field or opposite to it. If we let N+ and N respectively denote the number of dipoles oriented with 
and opposing the field, then the application of (2.2) yields 

(2.3) N + = Ce pB f kBT , iV_ = Ce ~ vElkBT . 


Since N = N+ + iV_ , it then follows that 


(2 - 4 ) N = 2Ccosh {Sr) 

The polarization for this configuration is given by 

p = pN + - pN~ = 2 pC sinh 



which yields the Ising spin relation 

(2.5) P = PN tanh(0) 


relating E and P. 

As detailed in [15], the anhysteretic polarization saturates to the value P s for increasing field inputs. 
Furthermore the relation (2.5) ignores the interaction with neighboring domains as well as electromechanical 


inputs due to applied stresses. 

( 2 . 6 ) 


The inclusion of these mechanisms 
Pan = Ps tanh 


(f 


yields the anhysteretic relation 
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where 


(2.7) 


E e = E + aP + P a 


denotes the effective field acting on the domain. The parameter a quantifies the degree of interdomain 
coupling while P a incorporates field contributions from an applied stress a. The parameter a quantifies a form 
of temperature-dependence due to the thermal energy [15]. For material characterization, the parameters 
P s ,a and a are estimated either from asymptotic relations or a least squares fit to data [16]. 

A second anhysteretic model is obtained under the assumption that the dipoles can orient uniformly in 
all directions. Integration and scaling for this case yields the Langevin model 


( 2 . 8 ) 


Pan — Ps 



As illustrated in [15], the Langevin model saturates less quickly than the Ising spin model since dipoles 
have more freedom concerning the directions in which they can orient. Both models have been employed to 
characterize the anhysteretic behavior of ferroelectric and piezo ceramic materials. 

The second component of the hysteresis model incorporates the energy required to translate and bend 
domain walls pinned at inclusions inherent to the material. As detailed in [15], this respectively yields an 
irreversible component Pi rr and reversible component P r ev to the polarization. The quantification of energy 
required to break pinning sites yields the differential equation 


(2.9) 


dPi ri 

~dE 


= 8 


Pan P ir 


k8 Ol (Pan Pirr) 


specifying the irreversible polarization. The parameter 8 = sign(dP) ensures that the energy required to 
break pinning sites always opposes changes in polarization. The physical observation that polarization 
changes after a reversal in field direction are reversible motivates the incorporation of the parameter 


j f 1 , {dE > 0 and P < P an } or {dE < 0 and P > P an } 
y 0 , otherwise 

Finally, the parameter A;, which quantifies the average energy required to reorient domains, is demonstrated 
in [16] to be asymptotically approximated by the coercive field E c in soft materials. 

The second component of the polarization is the reversible polarization which models the effects of 
domain wall bending. To first approximation, this is modeled by the relation 


(2.10) Prev — c (Pan Pirr) 

where c is a parameter which must be estimated for the specific application. 
The total polarization is then given by 


(2.11) P — Prev “b Pirr • 

To implement the model, the effective field for a given field and irreversible polarization level is computed 
using (2.7). This effective field value is then employed in either (2.6) or (2.8) to compute the corresponding 
anhysteretic polarization. The subsequent irreversible polarization is determined by numerically integrating 
(2.9) and the total polarization is specified by (2.11). We note that the prevailing polarization (2.11) is 
employed in (2.7) when modeling the contributions of neighboring dipoles on the effective field. 
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3. Time-Dependent Model. The model summarized in Section 2 is derived under the assumption 
of equilibrium conditions when employing the electrostatic potential energy relation (2.1) to derive the 
model for the anhysteretic polarization along with the energy required to reorient dipoles and hence break 
pinning sites. This assumption is valid under static or quasistatic operating conditions but omits rate- 
dependent mechanisms which are significant even at low frequencies. In this section, we derive a model for 
the anhysteretic polarization which incorporates this rate-dependence and reduces to the Ising spin model 
as the frequency is limited to zero. This quantifies a significant component of the rate-dependent behavior 
observed in Figure 1.1. 

3.1. Anhysteretic Polarization. To derive the anhysteretic model, we consider the material to be 
comprised of a lattice of cells with each cell having a dipole moment that is aligned either in the direction of 
the field or opposite to it. This is the same assumption made when deriving the Ising spin model in Section 2 
and is analogous to the regimes considered in [1, pages 69-71], [5, pages 104-106] and [19, pages 370-373]. 
The number of cells aligned in the direction of the applied field at time t is denoted by Ni(t) and N 2 (t) 
denotes the number of cells whose dipole moment is oriented in the opposite direction. 

The potential energy associated with the two equilibria is depicted in Figure 3.1. In the absence of an 
applied field, any interchange of dipole orientations is due to thermal fluctuations, whereas dipoles have a 
higher probability of overcoming the energy barrier and the number N\ ( t ) increases in the presence of an 
applied field E. To quantify this increase, we let W 21 denote the probability that one cell, considered for one 
second, switches orientation into the field direction and let W 12 denote the probability that a dipole switches 
in the opposite direction due to thermal excitation. The change in the number of cells having a specific 
orientation is then determined by the equations 


dNi 

dt 


= -^12^! -f W 2 \N 2 


dN 2 

dt 


— W12N1 — W21N2 . 


At equilibrium, the conditions ^ 
(3.1) 


^2- = 0 yield the requirement 

W12 _ N2 
W21 Ni 


In this case, the Boltzmann relation (2.2) holds and 


jVi = Ce pE/kBT 
N 2 = Ce~ pE/kBT 


where, from (2.4), the parameter C is specified by C(E,T) = 2 co&h(pE/k B T m y ^ comparison of (3.2) with 
(2.3) illustrates that the Ising spin distribution results when the time-dependent system considered here is 
allowed to reach equilibrium. Given the form of N\ and N 2 , the probabilities can be expressed as 

... = ±e--'“ T 

„„ = 

where r(E , T) is an arbitrary function of E and T. To determine the form of r(E, T), we note that when the 
system is not in equilibrium, and hence ^ 0, the probability that a cell switches orientation is nonzero. 
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Fig. 3.1. The double potential well associated with the two equilibria ; (a) No applied field , ( b ) An applied field E. 


Hence 

1 

W12 + W 2 1 — — 
T\ 


which yields the condition 

(3.3) t{E, T) = Ti (E, T) cosh (j^J ■ 

The parameter t\ is, in general, also a function of E and T. For the model comparison to PZT5A data 
presented in Section 4, we consider isothermal conditions and assume sufficiently small field dependence to 
justify considering T\ as constant. In general applications, however, the determination of mechanisms for 
determining the field and temperature dependence in T\ may improve the accuracy of the model. 

We now consider the solution of the resulting system 

— = -T e -pE/k B T N + ±_ e pE/k B T N2 

dt 2 T 2 T 

(3.4) 

dN 2 _ J -pE/k B T N _ ±_ e pE/k B T N 
dt 2 r 2 t 

under two sets of assumptions regarding the input field E. For a general E, the solution of (3.4) is 

Ni(t) = h+ jfc 2C -(‘/r)cosh(pE/fc B T) 


N 2 (t) = k l er 2pElkBT - k 2 er ( - t/T)cosh(pR/kBT) . 

To determine the integration constants k\ and k 2 , we consider the limiting behavior as t — ¥ 0 and t — > oo. 
As t — > oo, the distribution of dipoles limits to the equilibrium case (3.2) modeled by the Ising spin relation. 
This yields 

fci = Ce pE/kBT . 


The enforcement of equal initial dipole distributions, iVi(O) = A r 2 (0) , requires that 

k 2 =C sinh j 

The final distribution of dipoles at a given field level is then 


\k B Tj 


(3.5) 


Ni(t) = Ce pE/kBT + Csinh e_</Tl 

N 2 (t) = Ce~ pE/kBT - C sinh (j^f) e" t/n 
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where (3.3) was used to eliminate r(E,T). 

The polarization generated by this dipole configuration is specified through the relation 

P(t) = pN l (t)-pN 2 (t) 

(3.6) = 2pCsinh (iS ; )[ 1 ' e ’' / '‘] 

For the effective field E e = E + aP + P a introduced in (2.7), this yields the rate-dependent anhysteretic 
polarization model 

(3.7) Pan(t) = Ps tanh (^) [l - e~ i/n ] . 

A comparison of (3.7) with the Ising spin model (2.6) indicates that the latter is obtained as t — * oo. 

We consider now the response of the anhysteretic polarization to a periodic input field 

(3.8) E(t)=E 0 e iut . 


To accommodate the periodicity, we consider solutions of the form 

Ni (t) — No + v 0 e} ut 
N 2 (t) = N 0 - voe VjJt . 

The substitution of these expressions into (3.4) and consolidation of terms yields 


v 0 = 


( pE{t) \ 
1 + iojn ““““ V k B T ) 


N 0 


- tanh 


-iut 


where E(t) is of the form assumed in (3.8). The resulting polarization is then 


P(t) = pN^-p^it) 


pN 

1 + iiJTi 


tanh 



The isolation of the real component of the polarization, enforcement of saturation criteria, and the incorpo- 
ration of the effective field E e then yields the frequency-dependent anhysteretic relation 


(3.9) 


Pan{t) = 


1 +CJT1 


tanh 


(f 


We first note that as cj -> 0, this expression reduces to the Ising spin model (2.6). In this formulation, 
it can also be observed that the parameter t\ acts as a relaxation time. While T\ is considered constant 
for the example in Section 4, it is in general a function of both E and T. Finally, we note that for in- 
creasing frequencies, the expression incorporates the decrease in polarization observed in the data plotted in 
Figure 1.1. 
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3.2. Hysteresis Model. The expression (3.9) provides a frequency-dependent model for the anhys- 
teretic polarization. This relation is combined with the expression (2.9) for the irreversible polarization 
P irr and (2.10) for the reversible polarization P rev to obtain a frequency-dependent model for the total 
polarization. We note that the construction of the model in this manner neglects the quantification of rate- 
dependence in the energy required to bend domain walls or break pinning sites, thus rendering the model 
accurate only at low frequencies. The incorporation of these latter effects is under investigation and will be 
reported in a future work. An example illustrating the performance of the model at low frequencies (0.1 Hz 
and 1.0 Hz) is provided in the next section. 

4. Model Validation. The model developed in Section 3 quantifies both the hysteresis inherent to the 
E-P relation and the decrease in polarization which occurs when the frequency of the input field is increased 
from quasistatic to low frequency regimes. To illustrate its capabilities, we consider the characterization of 
polarization generated in a PZT5A wafer in response to a 2.5 MV/m input field at frequencies ranging from 
0.1 Hz to 1.0 Hz. A constant temperature was maintained to ensure isothermal conditions. 

The parameters a = 1.8 x 10 6 C/m, a = 3.6 x 10 6 Vm/C, k = 1.65 x 10 6 C/m 2 , c = 0.1, P s — 0.52 C/m 
and T\ — 0.16 were obtained through a least squares fit to the data collected at 0.1 Hz. Once obtained, these 
parameters were held fixed and variations in operating regimes were incorporated through the magnitude 
Eq and frequency lo of the input field. 

The model predictions at 0.1 Hz and 1.0 Hz are compared with the experimental data in Figure 3.2. It 
is observed that the model very accurately quantifies the data at 0.1 Hz which is the regime in which the 
parameters were estimated. The model also accurately predicts the decrease in polarization at 1 Hz, but does 
not yet incorporate the decrease in the coercive field which occurs as frequencj^ increases. The incorporation 
of rate-dependence in the energy required to bend and translate domain walls is under current investigation. 

5. Concluding Remarks. This paper addresses the quantification of certain rate-dependent mecha- 
nisms inherent to the hysteretic relation between the field applied to a piezoceramic material and the resulting 



Fig. 3.2. Model fit to 0.1 Hz and 1.0 Hz PZT5A data with the parameters a = 1.8 x 10 6 C/m, a - 3.6 x 10 6 Vm/C, 
k = 1.65 x 10 6 C/m 2 , c = 0.1, P s = 0.52 C/m 2 and n = 0.16. 
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polarization. The analysis presented here focuses primarily on the reduction in polarization which occurs as 
frequencies increase from quasistatic levels to low frequency regimes. This is modeled by determining the 
probability that dipoles achieve the energy required to overcome energy barriers and switch orientation when 
an external field is applied. The resulting model, which quantifies the anhysteretic polarization exhibited 
by the material, reduces to the Ising spin model when the driving frequency is reduced to quasistatic levels. 
This anhysteretic relation is then combined with expressions quantifying domain wall losses to provide a 
model which characterizes the hysteresis exhibited by the material. 

In its current formulation, the model is restricted to low frequency regimes since it does not yet incor- 
porate rate-dependent mechanisms for quantifying the energy required to reorient dipoles when translating 
domain walls. The extension of the model to include these mechanisms is necessary in order to increase the 
applicable frequency range of the model and is under current investigation. 
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